#ifndef _POSSIONOP_H_
#define _POSSIONOP_H_

#include "RegularGrid/RectDomain.h"
#include "Core/Tensor.h"
#include "Core/TensorSlice.h"

template <int Dim> class PossionOP;

template <int Dim>
class PossionOP{
public:
  using rVec = Vec<Real,2>;
  
  PossionOP(const RectDomain<Dim>& aDomain);

  void Laplacian(const Tensor<Real,Dim>& phi, Tensor<Real,Dim>& res);

  void relaxLaplacian(const Tensor<Real,Dim>& phi, const Tensor<Real,Dim>&
                      rhs, Tensor<Real,Dim>& res, Real weight);
  
  void computeResidul(const Tensor<Real,Dim>& phi, const Tensor<Real,Dim>&
                      rhs, Tensor<Real,Dim>& res);
  
  Real computeNorm(const Tensor<Real,Dim>& phi, int p = 0);
protected:
  RectDomain<Dim> Domain;
};

template <int Dim>
PossionOP<Dim>::PossionOP(const RectDomain<Dim>& aDomain){
  assert(aDomain.getCentering() == NodeCentered);
  assert(aDomain.getNumGhost() == 1);
  Domain = aDomain;
}

template <>
void PossionOP<2>::Laplacian(const Tensor<Real,2>& phi,
                             Tensor<Real,2>& res){
  const Box<2>& bx = Domain;
  const rVec& dx = Domain.spacing();
  assert(res.box().contain(bx));
  loop_box_2(bx,i,j){
    res(i,j) = (-2*phi(i,j)+phi(i-1,j)+phi(i+1,j))/(dx[0]*dx[0]) +
      (-2*phi(i,j)+phi(i,j-1)+phi(i,j+1))/(dx[1]*dx[1]);
  }
}

template <>
void PossionOP<1>::Laplacian(const Tensor<Real,1>& phi,
                             Tensor<Real,1>& res){
  const Box<1>& bx = Domain;
  const Vec<Real,1>& dx = Domain.spacing();
  assert(res.box().contain(bx));
  loop_box_1(bx,i){
    res(i) = (-2*phi(i)+phi(i-1)+phi(i+1))/(dx[0]*dx[0]);
  }
}

template <>
void PossionOP<2>::relaxLaplacian(const Tensor<Real,2>& phi, const Tensor<Real,2>&
                                  rhs, Tensor<Real,2>& res, Real weight){
  const rVec& dx = Domain.spacing();
  const Box<2>& bx = Domain;
  assert(res.box().contain(bx));
  const Real b = 2.0/(dx[0]*dx[0])+2.0/(dx[1]*dx[1]);
  loop_box_2(bx,i,j){
    Real tmp = (phi(i-1,j)+phi(i+1,j))/(dx[0]*dx[0])+(phi(i,j-1)+phi(i,j+1))/(dx[1]*dx[1]);
    Real r = (tmp - rhs(i,j))/b;
    res(i,j) = r * weight + phi(i,j) * (1 - weight);
  }
}


template <>
void PossionOP<1>::relaxLaplacian(const Tensor<Real,1>& phi, const Tensor<Real,1>&
                                  rhs, Tensor<Real,1>& res, Real weight){
  const Vec<Real,1>& dx = Domain.spacing();
  const Box<1>& bx = Domain;
  assert(res.box().contain(bx));
  const Real b = 2.0/(dx[0]*dx[0]);
  loop_box_1(bx,i){
    Real tmp = (phi(i-1)+phi(i+1))/(dx[0]*dx[0]);
    Real r = (tmp - rhs(i))/b;
    res(i) = r * weight + phi(i) * (1 - weight);
  }
}

template <int Dim>
void PossionOP<Dim>::computeResidul(const Tensor<Real,Dim>& phi, const Tensor<Real,Dim>&
                                  rhs, Tensor<Real,Dim>& res){
  Laplacian(phi,res);
  res = rhs - res;
}

template <int Dim>
Real PossionOP<Dim>::computeNorm(const Tensor<Real,Dim>& phi, int p){
  const Box<Dim>& bx = Domain;
  if(p == 0)
    return norm(phi.slice(bx), 0);
  else if(p == 1)
    return norm(phi.slice(bx), 1) * prod(Domain.spacing());
  else if(p == 2)
    return norm(phi.slice(bx), 2) * sqrt(prod(Domain.spacing()));
  return 0.0;
}


#endif //_POSSIONOP_H_
